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Abstract 

In two previous works [arXiv:1009. 4363, arXiv:1107. 0668], we studied the time evolution 
of a system of real scalar fields with quartic coupling which shares important features with 
the Color Glass Condensate description of heavy ion collisions. Our primary objective was 
to understand how such a system, when initialized with a non-perturbatively large classical 
field configuration, reaches thermal equilibrium. An essential goal of these works was to 
highlight the role played by the quantum fluctuations. 

However, these studies considered only a system confined within a box of fixed volume. In 
the present paper, we extend this work to a system that expands in the longitudinal direction 
thereby more closely mimicking a heavy ion collision. We conclude that the microscopic 
processes that drive the system towards equilibrium are able to keep up with the expansion 
of the system; the pressure tensor becomes isotropic despite the anisotropic expansion. 

1 Introduction 

The issue of thermalization of the quark-gluon matter produced in heavy ion collisions is one 
of the most challenging problems in this field. On the one hand, this matter appears to behave 
like a nearly perfect fluid, as is suggested by the comparison between flow measurements at 
RHIC [1-4] and hydrodynamical models with low shear viscosity [5]. Moreover, this comparison 
seems to call for a very early onset (< 1 fm/c) of the hydrodynamical behavior. Such a good 
description by hydrodynamics is usually interpreted as an indication that the system is fairly close 
to local thermal equilibrium. On the other hand, trying to justify this fast equilibration from first 
principles has proven to be notoriously difficult and no definitive conclusion has been reached 
thus far. Early estimates of the thermalization time based on a standard kinetic description (with 

2 — !• 2 and 243 processes), as in the bottom- up scenario [6], lead to equilibration times that 
are much larger than what hydrodynamics seems to require. In addition to kinetic scattering 
processes, quarks and gluons with an anisotropic particle distribution are subject to Weibel 
instabilities [7, 8] that are well known in QED plasmas. These potentially play a significant role 



in restoring isotropy and in thermalizing the system. A formal description of Weibel instabilities 
can be couched in the framework of Hard Thermal Loop effective field theory where the hard 
particles coexist with soft fields [9-13]. This approach has been applied to study the evolution of 
such anisotropic systems, and the role of instabilities in their relaxation towards equilibrium, in a 
series of mostly numerical works [14-24] . For recent analytical parametric estimates of the effect 
of Weibel instabilities in the expanding systems produced in heavy ion collisions, see refs. [25, 26]. 

Another feature of heavy ion collisions that complicates our understanding of thermalization 
is the fact that, at such early times, the collision is more appropriately described in terms of fields 
rather than on-shell particles 1 . Such a description arises naturally in the Color Glass Condensate 
(CGC) effective theory. In this approach [27-33], the nuclei are described in terms of highly 
boosted color sources coupled to soft gauge fields. At the very high energies reached at RHIC and 
LHC, the density of color sources becomes so large that their currents are inversely proportional 
to the (assumed to be small) coupling constant. In this framework, one can systematically arrange 
the various contributions to a computation of inclusive quantities in powers of the strong fine 
structure constant a s with a given order in a s corresponding to a fixed loop order. The leading 
order (LO) approximation can be expressed entirely in terms of the classical solution of the 
Yang-Mills equations, the next-to-leading (NLO) order is composed of the 1-loop graphs in the 
presence of a classical background field, and so on. A very satisfying aspect of this description is 
that the color source distribution of the projectiles appears only via universal weight functionals 
of these distributions. Factorization theorems have been proven that demonstrate (to leading 
accuracy in powers of a s Y , where Y is the rapidity) that these weight functionals arc the same 
for all inclusive observables and for different reactions involving identical projectiles [34-36]. The 
evolution of these weight functionals with rapidity is described by the JIMWLK equation[37-44]. 

In the CGC framework, the Weibel instabilities manifest themselves in the form of unstable 
solutions of the classical Yang-Mills equations [45-53], that lead to secular divergences in the 
NLO 2 correction to quantities such as the pressure: these corrections diverge when the time 
goes to infinity. In ref. [57], we sketched an improvement that resums, at each order of the 
expansion in a s , the terms that have the fastest growth in time. This resummation scheme, 
which amounts in practice to averaging classical solutions of the field equations of motion over 
a Gaussian ensemble of initial conditions 3 , was justified more extensively in [65, 66]. Moreover, 
the precise form of the initial Gaussian ensemble of gauge fields was also derived in this work. 

Since the numerical implementation of this resummation in QCD is computationally very 
demanding, we first considered a much simpler field theory in [67, 68], a scalar field theory with 
a quartic <f> A coupling. Similar issues arise in this simpler context and provide one with a feeling 
for the effect of instabilities on the evolution of the system towards equilibrium. The scalar 
theory, like QCD, is scale invariant at the classical level in 3 + 1 dimensions 4 . More importantly, 
this theory is known to have unstable classical solutions because of parametric resonance (see 
refs. [69, 70]). Although the microscopic nature of these instabilities is very different from the 
Weibel instabilities encountered in Yang-Mills theory, they also lead to secular divergences in 
the pressure at NLO. Therefore, as in the latter, a resummation of the divergent contributions 
at each order is essential to obtain finite results at late times. 

In refs. [67, 68], a resummation identical to the one advocated in [57, 65, 66] for QCD 

-'-Soft modes, as we know from the uncertainty principle, need a rather long time to go on-shell. 

2 At LO, the pressure is finite at all times, but the longitudinal pressure is negative [54-56]. 

3 To the best of our knowledge, this scheme was first used in [58-60], in the context of preheating in post- 
inflationary cosmology. See also ref. [61]. In the context of non-abelian plasma instabilities, it was used already 
in [62], A similar approach has also been applied to cold atom physics, in problems related to Bose-Einstein 
condensation (see for instance refs. [63, 64]). 

4 This property is not mandatory, but it has simple implications for the equation of state: e = 3P, up to small 
corrections. 
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was applied to this scalar theory and implemented in a numerical lattice computation. These 
simulations showed that the pressure quickly relaxes to the value required by the equilibrium 
equation of state (P = e/3). On longer timescales, the particle distribution evolves to a classical 
equilibrium distribution. Moreover we found that when the system is initialized with an excess of 
particle number (compared to the equilibrium value at a given energy), then a transient chemical 
potential appears, and even a Bose-Einstein condensate if the particle number excess is too large 
(see also refs. [71-73]). 

In the two papers already mentioned exploring non-equilibrium scalar dynamics, only systems 
contained in a fixed volume box were considered. The isotropization of the pressure tensor was 
therefore not relevant. To address this issue, we shall extend here our study of scalar </> 4 field 
theory to a system that expands longitudinally. In this case, one expects naively that the 
longitudinal expansion produces a red-shifting of the longitudinal momenta thereby leading to a 
small longitudinal pressure that decreases faster than the transverse one. If this conclusion was 
right, hydrodynamics would not be applicable to the description of such an expanding system. 
The main question we want to address with this extension is whether the instabilities can keep 
the longitudinal pressure close to the transverse one, despite the longitudinal expansion. A study 
of these instabilities for an expanding system has been presented in ref . [74] , but the main focus 
of this work was not on isotropization. 

We will mimic closely, in this scalar theory framework, the Color Glass Condensate description 
of high energy heavy ion collisions by assuming that the classical background field is independent 
of rapidity. The fluctuations that are superposed on this classical field at the initial time are the 
scalar analogues of the spectrum derived for QCD in [65, 66]. The derivation of this spectrum 
in the much simpler scalar case is done in section 2. In the section 3, we discuss the lattice 
discretization of the model. The time evolution of the particle distribution is studied in the 
section 4. This discussion also shows why it is important to have a very small lattice spacing in 
the rapidity direction. The time evolution of the transverse and longitudinal pressures is studied 
in the section 5. We conclude that the instabilities are indeed able to make the longitudinal 
pressure nearly equal to the transverse one. In the section 6, we compare the results obtained 
here in classical statistical field theory with expectations from viscous hydrodynamics. Section 7 
is devoted to conclusions and to an outlook into the application of this approach to QCD. Some 
technical notes and complements are relegated to three appendices. 



2 Spectrum of fluctuations in an expanding system 

Following refs. [67, 68], we consider a massless real scalar field theory with a quartic coupling, 
whose Lagrangian density is given by 

£(0) = ^0) 2 -f^ 4 . (i) 

The kinematics of a high energy nuclear collision is mirrored by picking a preferred spatial 
direction - the z direction in this paper - to be the collision axis. For collisions in the high 
energy limit, the problem is invariant under finite Lorentz boosts in this direction. A natural 
system of coordinates in which this invariance takes a simple form are the proper time r and the 

5 These two features can only be transient because the particle number is not conserved in this theory and 
eventually inelastic processes will eliminate any excess in the particle number. 
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rapidity r\ coordinates that are related to the usual Cartesian coordinates by the transformations 



with the transverse coordinates x and y unchanged. In this coordinate system, the d'Alembertian 
reads 

□ = d 2 T + ~d T ~ \tf - Vi . (3) 

r t 

Longitudinal boosts preserve t and x±, and shift the rapidity by a constant. Thus, the invariance 
under longitudinal boosts implies that physical quantities are independent of the rapidity r/. 



2.1 Classical background field 

In the CGC description of heavy ion collisions, the leading order approximation is given by the 
retarded classical solution of the Yang-Mills equations. In subsequent orders, this classical field 
plays the role of a background field whereby all the propagators and vertices are dressed. This 
classical background field is boost invariant. Thus, in our scalar analog, we will require the 
background field to be independent of r\. The classical equation of motion for such a field ip is 

<p+ 1 <p-ViLp + V'(Lp)=0, (4) 

T 

where the dot denotes a derivative with respect to the proper time r. Note that we only consider 
here the evolution of the system after the collision (for r > 0) . Because the sources that describe 
the colliding projectiles have support only on the axis z = ±t, the field evolution is not driven 
by these sources; instead, they provide initial conditions for the field at r = + . 

To ascertain the nature of the initial conditions at r = + , let us neglect for the time being 
the interaction term V'(ip) and decompose the x± dependence of the field in Fourier modes as 

<P(r,x ± ) = J g^ fcx (r) e *— . (5) 

The Fourier coefficients obey the equation 

<fk ± + -<j>k± + fci^fcj_ = . (6) 

T 

This is a Bessel equation of index n = and its general solution is of the form 

<Pk ± (r) = a k± Jo(k±r) + b k± Y {k±T) , (7) 

where a k± and b k± are constants. Only the term that contains Jo can smoothly approach the 
limit t — > + . (If we keep the Y term, the field would diverge as ln(l/r) and its derivative as 
1/r in this limit.) 

Therefore the solutions that are regular near the origin are of the form, 

<P{r,x ± ) =J^2 etk± ' X± «*x Mk±r) , (8) 
where the a k± coefficients are just the Fourier transform of the field at r = 

<p{0 + ,x ± )= J ^j k± - m± Ofcx- (9) 
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Note that its derivative tends to zero at the origin, 

lim (f(T,x±) = . (10) 

T->0+ 

The interaction term V'(ip) that we neglected in this discussion introduces some mixing 
among the Fourier modes and modifies the evolution of the solution. However this modification 
is negligible at short times because the evolution of the field at very short times is dominated by 
the expansion rate which is much larger than the interaction rate. Therefore the solution that 
we find by neglecting the interaction term remains valid in the immediate vicinity of the origin 6 . 



2.2 Spectrum of fluctuations 

We want now to perform for this scalar model the resummation described in ref. [65, 66] for the 
QCD case. This resummation amounts to allowing the initial condition for the classical field to 
fluctuate and to average observables over these fluctuations. Following [66], we can construct the 
Gaussian ensemble of initial conditions by superposing on the background field ip (independent 
of rj) a fluctuating n dependent term, with the result expressed as 



<f> = <P+ 



+ c a 

K K K K 



(11) 



Here the a K are mode functions propagating on top of the classical background field tp that obey 
the equation of motion 



1 1 

-a K ; 

T T' 



" ■ - — d n a K - V l a K + V"(<p)a K = . (12) 



This is a linear equation and the set of its solutions is a vector space. The symbol K denotes 
collectively all the labels that are necessary in order to index the basis for this vector space and 
dfx x is the measure necessary to integrate over this space. In eq. (11) the a K 's must be positive 
frequency solutions and they must form an orthonormal basis with respect to the inner product 

(a\b) = it J dr)d 2 x x (a*b-a*b) . (13) 

This means that one must have 

(OkK) = S KL . (14) 

where the (5-function on the space of the indices K is normalized to ensure 



J dfx x 5 KL = 1. (15) 



The coefficients c K are Gaussian random complex numbers whose variance is defined by the 
equations, 

i c K c L ) = («)=0 

(c K c* L ) = \S KL . (16) 

We emphasize that the normalization condition of the fluctuations a K depends on how we define 
the integration measure so that at the end of the day the ensemble of <p's defined by eq. (11) 
does not depend on this arbitrary choice. 

6 How far in r this approximation is valid depends on the magnitude of the field and of the coupling constant 
g, both of which determine the strength of the neglected non-linear terms. 
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2.2.1 Mode functions 

Let us now solve eq. (12) which determines the mode functions. We wish only to obtain its 
solutions at small times r, therefore we shall neglect the time dependence' of the background 
held ip and replace it by its limit at r — > + , 

fo(x±) = lim <p(r,x±) . (17) 

T-S-0+ 

Since ipo does not depend on the rapidity rj, we can look for solutions that have a well defined 
wave-number in the variable r\ of the form 

a„(T,r),x ± ) = e ivri b v (T,x ± ). (18) 

The function b v now obeys 

1 v 2 

K + -K + -zK - V\b u + V"(<po)K = . (19) 

We note that the x± dependence of the solution is controlled by the operator — + V"(ipo). 
This operator is real and symmetric and can therefore be diagonalized over a basis of orthogonal 
functions Xk( x ±), 

J d 2 x x xl(x±)xi(xi_) = 5 kl . (20) 

The ^-function S k i is dehned so that its integral over the measure d/i k in fc-space is 

d/-ifc S k i = 1 ■ (21) 

The next step is to decompose the x± dependence of b v on the basis provided by the Xk( x ±Y s - 
This amounts to looking for fluctuations that are of the form 

a vk (T, 77, x±) = e lvn Xk(x±) &„fc(r) . (22) 

The remaining function b vk {r) depends only on time and must obey the following equation 

1 • v 2 

Kk + -b v k + ~^b vk + X k b vk = . (23) 

T T 

II JUIJ Wilt ,1^1 ' I H ' I f I J M'1UMII!I lllll'fll. \ 1 1 1 1.1 I II II fl I K I 1 I III Mil' llillilMVJ IIII.H I HHl^ 

H^l\\ k T) and H^ 2 \x k r). The latter is the positive frequency solution, as can be seen from its 



This is a Bessel equation whose general solution is a linear combination of the Hankel functions 

it1 

asymptotic behavior at large time 9 . Thus the a K 's in eq. (11) can be chosen as 



<Wt,?7,xx) = OL vk e ivr > X k(x±) h£\\ht) . (24) 



7 Since we saw in the previous section that its time derivative should vanish when r — > 0+ . 

8 One could also use the Bessel functions Jiv(\k T ) an d Yiv(^k T ) as the two independent solutions, but they 
are less convenient in the problem at hand because they are both mixtures of positive and negative frequency 
components. 

9 For large arguments, one has 

(2), v -z( T -inv/2-n/4) 



G 



We need only to determine the constant prefactor a vk in order to ensure that these solutions 
are properly normalized. (They are, by construction, already mutually orthogonal.) To obtain 
this, we can compute the inner product 



= it 2tt5(v - fx)6 kl |a„ fc | 2 H^*(X k T) Or Hg\\ kT ) 

v ' 

— 4ie — 77 v / ttt 

= 86(v-n)5 h i K fc | 2 e""" , (25) 

where the identity used to go from the second to the third line is the Wronskian between — 

[H^y and . If we choose the integration measure over the indices v, k to be dji K = d/Zfc, 
then this inner product should be 2ir8{v — (x)S k i, which can be satisfied by taking 

a vk = . (26) 



2.2.2 Final form of the spectrum 

Let us now summarize the results of this section. We showed that the fluctuating fields defined 
in eq. (11) should have the explicit form 

0(t, r,, sl) = p(r, x ± ) + ^-J^ d ^ e™' 2 c vk e*"> X k(xi_) H^(X k r) + c.c. , (27) 

where the random coefficients in this linear superposition have the variance 

(cvkc^i) = , (c„ fe c*j) = u8{v - n)Ski ■ (28) 

Thus far, we haven't been very explicit about the nature of the index k that labels the eigen- 
functions of the operator — + V"(ipo). This label is a 2-dimensional continuous index that 
plays the same role as the transverse momentum in the free case 10 . 

The fields defined by eq. (27) do not seem to have a well defined limit when r — » + , because 
the Hankel functions oscillate like T ±tu as r approaches and their derivatives diverge. This 
forces us to start the time evolution at a small but finite time To, rather than at exactly t = + . 
In this respect the fields with rapidity dependent fluctuations superposed on the background field 
are qualitatively different from the background field, for which starting from r = + poses no 
difficulty. However, the finite To that we are forced to introduce as a consequence, fortuitously 
has no effect on the physical results. Indeed since eq. (27) is a solution of the equations of motion 
at small r, one can prove that physical quantities computed at some later time by averaging over 
these 0's will be independent of the choice of To provided it is chosen to be small enough. This 
will be checked in the section 3.3. 

10 In the free case, the corresponding measure and delta function would thus be defined as 

d 2 fcx 

d A'fe = , n x9 > & kl = (27r) 2 <5(fc x - l ± ) . 
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3 Lattice implementation 



An explicit determination of the fields in eq. (27) requires that we solve the eigenvalue equation in 
eq. (20). Because this expression depends nonlinearly on the initial background field, the solution 
has to be determined numerically. In this section, we will outline the numerical procedure followed 
in the solution. 

3.1 Discretization 




Figure 1: Details of the lattice discretization used in our numerical implementation. 
In view of a numerical implementation, we discretize the problem as follows: 

i. The transverse plane x± becomes a L x L lattice, and the lattice sites are labeled by a pair 
of integers (that range from to L — 1). We use periodic boundary conditions. The 
indices i, j are thus defined modulo L. We follow the convention to set the transverse lattice 
spacing to unity, a± — 1. All the other dimensionful quantities of the problem are then 
expressed in units of the appropriate power of the transverse lattice spacing. For example, 
the values of the time quoted later on in the paper should be understood as values in units 
of t/ ax- 
ii. We consider a unit slice of rapidity, discretized in N intervals, labeled by an integer n that 
ranges from to N —1. We also adopt periodic boundary conditions here. This is justified 
for high energy collisions because we expect observables to be invariant by translation in r\ 
for the r\ range of interest for thermalization. The index n is thus defined modulo N. The 
lattice spacing in the rapidity direction is dimensionless and its value is h = 1 /N. 

iii. We discretize the second derivatives in 77 and x± by symmetric finite differences as 

d%(f> -> {4>n+lij + <f>n-Uj - 2<f>nij) 

(29) 
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iv. Time remains a continuous variable 11 . 



On the lattice, the eigenvalue problem for finding the Xfc' s is transformed into the problem 
of finding the eigenvalues and eigenvectors of a real symmetric matrix, since 



[-Vi + V"(<p )] X -> (Dx)ij = ^Xij ~ Xi+ij ~ Xi-ij ~ Xij+i - Xij-i + V-jXij ■ 
Thus, the L 2 x L 2 matrix T)^ ^ that we must diagonalize has the following expression 

Dij t ki = (4 + V-j)5 ik Sj[ - (5 i+ ik + 6i-ik)6ji - 5 ik {5j+u + . 
This diagonalization will provide us a set of L 2 eigenvectors such that 

E D *iM Xw } 

kl 

Vy (p) *y (9) = L 2 5 



\ 2 

A(p)Xij 



(30) 



(31) 



(32) 



Once this is done, the fields of eq. (27) can be written as 



2NL 2 h 



h E eW2 E [«* ^ *S? flSw +c.c. 

p=0 



v=0 



where v is an eigenvalue of the discretized operator d 2 

, 2 



h 2 



1 — cos 



/2tto\ 

{ir) 



and where the c vp are random Gaussian numbers with the variance 



(c c* )-^6 S 

y-n-p'-raql ~ r> VnmVpq ■ 



(33) 



(34) 



(35) 



3.2 Special case of a uniform background field 

If the background field tp(r, x±) is independent of x±, the diagonalization of the matrix Dij^i 
is trivial because VZ- is just a constant mass term m 2 = g 2 y ^° 2 ^ , independent of position in the 
transverse plane. The eigenfunctions are thus plane waves labeled by two integers p, q between 
and L — 1, 

Xij = e L , (36) 

and the corresponding eigenvalues are 



\2 



m 2 + 2 



2 — cos ( 



(?) 



(?) 



(37) 



11 Time steps in the numerical implementation can be made as small (and even changed dynamically) as needed 
to ensure the desired accuracy. 
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3.3 Independence with respect to the initial time 



As noted previously, one cannot take the time to zero in equation (27) due to the oscillatory 
behavior of the Hankel functions. Thus one must choose a small initial time To, compute the 
initial fields at this time, and proceed with the time evolution from there. However, since this 
initial time is arbitrary, observables should not depend on it. One could view this question 
within the framework of the renormalization group: observables can be made independent of tq 
by making the ensemble of initial fields depend on To in such a way as to compensate for the fact 
that the time evolution starts at t . In fact, the time dependence of eq. (27) is precisely what is 
needed to achieve this. 

In figure 2, we show the time evolution of the energy density e and transverse pressure Pt 
(see eqs. (48) for explicit expressions in terms of the scalar field) of the system for two different 
initial times, tq = 10 -2 and Tq = 10 _1 . One can see that, for the times that are common to 
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Figure 2: Time evolution of some components of the energy-momentum tensor, for two different 
initial times To = 0.01 (solid lines) and 0.1 (points). This computation was done on a 20 x 20 x 20 
lattice with 256 field configurations. All the plots presented in this paper (with the exception of 
fig. 17 in appendix B) correspond to coupling constant g = 4. 

the two evolutions, these physical quantities do not depend on the time at which the system is 
initialized. One should add a word of caution here. This statement is true only as long as the 
initial time To is small compared to the physical time scales in the problem (for example, the 
period of the oscillations of the pressure in the fig. 2). Note also that the classical background 
ip itself in eq. (27) evolves from To to Tq, although this is a small effect if both initial times are 
small. 

3.4 Unstable modes 

A crucial aspect of dynamics in this theory is the presence of unstable modes that are amplified 
by parametric resonance (see also ref. [74] for a discussion of this question). We saw in [67] that 
the resonance band for a constant volume is fixed and that its boundaries are proportional to 
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the amplitude of the background classical field. The situation becomes more complicated when 
the system is free to expand in the z direction. 

For simplicity, let us assume that the background field tp is spatially homogeneous 12 . The 
linearized equation of motion for a small fluctuation of wave numbers v and k± propagating over 
this background reads 




When the background field is not present (p = 0), the solutions of this equation are the Bessel 
functions Jiv(k±r) and Y iu (k±r). At late times, they oscillate with an amplitude that decreases 
like r" 1 / 2 . When tp ^ 0, one can check (numerically) at any given time whether the fluctuation 
of wavenumbers v, k± has the expected magnitude, namely, if they were decreasing at the pace 
of the above Bessel functions, or if they are larger than expected. If the latter is the case, one 
interprets this mode to have been amplified by the resonance in the course of its evolution. 
In fig. 3, we plot modes that are unexpectedly large at the time r for three values of r. The 
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Figure 3: Modes that have been amplified by the instability at various times in the evolution. 

number of modes that have been affected by the resonance increases with time. Moreover, the 
boundary of this region is well reproduced by a line of constant k\ + v 2 /t 2 . This fact is easy 
to understand semi-quantitatively from what we know about the resonance in the fixed volume 
case. There the resonant modes form a narrow band 

R- < k\ + k 2 z < R+ , (39) 

where R± are two numbers proportional to the amplitude squared of the background field. 
The resonance occurs due to a special tuning between the "effective mass" k 2 of a mode and 
its coupling to the oscillating background field. In the expanding case, it is the combination 
k\ + v 2 /t 2 that plays the role of a mass term for the mode. For large times such that the 
variations of 1/r 2 are slow compared to the oscillations of the background field we can justly 

12 This choice is the one for which the structure of the resonance band is the simplest. See for instance [69]. 
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transpose to this case the resonance condition obtained in the non-expanding case. We then get 
the resonance condition, 

v 2 

R- <k\-\ — ^ < R+ . (40) 

These inequalities define a domain that has the shape of a narrow corona in the fcj_ plane. 
The main effect of the expansion is that this corona expands with time, starting very close to 
the v = axis and progressively moving towards larger values of v. Therefore the modes that 
are unstable are not the same at all times. Moreover, a given v, k± mode can only be resonant 
during a finite amount of time while the resonance band is sweeping through it. During this 
time, it is amplified by some possibly large but always finite factor. 

From this discussion it is easy to understand fig. 3. The shaded areas represent all the modes 
that have been resonant at some point between the initial time and the time r of interest. This 
domain is the addition of the locations of the resonance band for all the times up to r. In this 
interpretation, the boundary of the domain (of equation k\ + v 2 /r 2 = const) corresponds to 
the modes that are just starting to become resonant at the time r. Note that the modes that 
become resonant at late times are not amplified as strongly as those that resonate early: indeed, 
as known already for a fixed box, the Lyapunov exponent is proportional to the amplitude of the 
background field-which in this case decreases with time. 



4 Occupation number 

In the next two sections, we will discuss results from our lattice simulations. In this section, we 
will discuss results for the occupation number. In the following section, results for components 
of the stress-energy tensor will be presented. 

4.1 Definition 

We begin with the Fourier decomposition of a free scalar field operator 13 in the system of coor- 
dinates (r, ?7, x±) 

?(r,f7,«±) = ^J e^H^{k^r)a vk± e** e*— + h.c. . (41) 

One may check immediately that this expression is a solution of the free field equation of motion. 
This equation can be inverted to obtain creation and annihilation operators expressed in terms 
of the field as 

a vk± = lT e™l 2 ^ J d 2 x ± d V e ~ ik ^ e** (k x r) d T 4(r,n,x ± ) . (42) 

The normalization in these formulas has been chosen to satisfy 

«»k ± ,al l± ] = {2n) 3 5(u - v')8{k x - Z x ) , (43) 

given the canonical commutation relations for the field operator. When we set v = v 1 and 
k± = l±, the delta functions on the right hand side become S±L V , where S± is the transverse 
area of the system and L n the length of the rapidity interval under consideration. 



3 The seemingly peculiar normalization factors in this formula will become clear after eq. (43). 
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As explained in ref. [68] (sections 4.1 and 4.2), it is straightforward to obtain the expectation 



value of the symmetric combination <va 



in the resummation scheme we are using here. 



For creation and annihilation operators normalized as in eq. 
occupation number as 



(43), it is natural to define the 



1 + 2f uk± = 



1 



S±L V 



(44) 



From eq. (42), the occupation number computed in our resummation scheme can be expressed 



as 



Uk ± (t) 



1 TTT^e 



2 AS^Lr, 



d 2 x ± d V e - ivr >e- ik ^ Hg>*(k±T) d T <f>(r,r),x± 



(45) 



Here </>(t,t],x±) is a classical solution of the field equation of motion and the brackets denote 
the average over the Gaussian ensemble of initial conditions. 

A crucial aspect of eq. (45) is the prefactor r 2 e Triy in front of the integral that ensures that 
each mode is weighted correctly at all times despite the dilution due to the longitudinal expansion 
of the system. One can also check that this formula gives a vanishing occupation number in the 
case of pure vacuum fluctuations (i.e. with <fi replaced by eq. (27) with ip = 0). Moreover, we 
have checked that, when evolved in time with the interacting equation of motion, this vacuum 
ensemble of fields evolves in such a way as to give a vanishing occupation number over all the 
time range of interest 14 . 



4.2 Time evolution 

In fig. 4, we show the occupation number at two times, computed on a 40 x 40 x 80 lattice, for 
an initial classical background field <fo(x±) at t = + . As an example, we take a background 
field that initially contains only a single mode fc_L, 

fo(x±) = (f x cos(fc_L • x±) . (46) 

This is arguably a simplistic toy model of the classical background field. In the Color Glass 
Condensate description of a heavy ion collision, the classical background field would contain a 
range of non-zero modes up to the saturation momentum. The coupling constant used in this 
computation, and in all the plots shown in this paper, is g — 4. The value of the background field 
has been taken to ipo = 15, and k± = 0.77. At the time r = 0.1, a harmonic of the initial k± 
mode has already started to grow in the distribution (see the second peak in the left plot of the 
figure 4) , but its amplitude is comparatively very small (note that the vertical axis is logarithmic) . 
On the other hand, the distribution in v is still confined at v = because the fluctuating part 
of the fields is still very small compared to the boost invariant classical background. The plot 
on the right of the figure 4, that represents the occupation number at the time r = 200, shows 
that a continuum of fcj^ and v modes eventually get populated after some evolution. 

This plot also illustrates the main difficulty in simulating an expanding system. Over time the 
distribution extends in the v direction and eventually reaches the lattice cutoff in this variable. 
The origin of this behavior can be easily understood. For simplicity, let us assume that the system 
has reached local thermal equilibrium and expands following ideal Bjorken hydrodynamics. Then 
its energy density decreases like e ~ r~ 4 / 3 and its temperature decreases like T ~ t" 1 / 3 . The 
temperature is roughly the extent of the support of the occupation number in the variables k± 

14 If this evolution is carried to extremely large times, we may expect that these vacuum fluctuations will 
eventually "thermalize" [75]. 
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T = 1 [40 x 4(1 x mi 



T = 200 [40x40x81) 




Figure 4: Occupation number at times r = 10 1 and r = 200, on a 40 x 40 x 80 lattice. 



and k z . Therefore in these variables the support of the occupation number tends to shrink slowly 
with time. However at late times the variables k z and v are related by 

k z ~ - , (47) 

T 

and therefore a shrinkage as r -1 / 3 of the momentum k z corresponds to an increase with t 2 / 3 of 
the variable v. Therefore if the system thermalizes 15 or approaches thermal equilibrium, there 
will always be a time where the support of f u k± hits the lattice cutoff in the variable v no matter 
how high this cutoff is. 

This is a serious issue for our simulation because physical quantities computed after this 
time will be contaminated by lattice artifacts. This is the case in particular for components of 
the pressure tensor as they are sensitive to the hardest populated modes in the system. When 
the occupation number reaches the upper v limit, the distribution in k z is artificially cut-off, 
which leads to a reduction of the longitudinal pressure. If we were to pursue the time evolution 
much beyond this time, the system would effectively become 2-dimensional and we would get 
the incorrect result P L <C P±. 

A computationally expensive but otherwise straightforward way to avoid this problem is 
to have a finer discretization in the rapidity direction, dividing the unit rapidity interval we 
are considering in a large number N of lattice spacings. As N increases, we can push the 
time evolution to later times. In the table 1, we give an estimate (obtained by computing the 
occupation number at various times and observing when its support first reaches v max ) of the 
maximal time, for a given N, when lattice artifacts are negligible. In the rest of the paper, 



N 


Tmax 


80 


120 


160 


220 


320 


300 



Table 1: Maximal time as a function of the number of lattice spacings in the r\ direction, 
we shall present results obtained on a 40 x 40 x 320 lattice, which allows us to go safely up to 

15 In contrast, if the system is free-streaming, kj_ ~ const and k z ~ r _1 . This corresponds to a constant support 
both in k± and in v. 
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r = 300. The occupation numbers at the times r = 10, 50, 100 and 300 are shown in the figure 
5. 

T=10 [40 x 40 x 320] T = 50 [40x40x320] 




Figure 5: Occupation number at times r = 10, 50, 100 and 300, on a 40 x 40 x 320 lattice. 



5 Energy-momentum tensor 
5.1 Definition 

We are now ready to explore the evolution of the energy-momentum tensor. Besides the relax- 
ation of the system towards an equation of state, the longitudinal expansion means that there is 
no guarantee that the transverse and longitudinal pressures will be equal. However, for hydrody- 
namics to be applicable, the anisotropy of the pressure tensor should not be large. It is therefore 
of the utmost importance to determine whether the longitudinal and transverse pressures ever 
approach each other for the expanding system. 

In the (r, rj,x±) system of coordinates, the diagonal components of the energy-momentum 
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tensor are given by 



T 



j*yy — 



'XX 



£ 



-(0 2 + (v ± 0) 2 + r- 2 (a t ,0) 2 )+y(0) 

i(0 2 -(a^) 2 + (a^) 2 -r- 2 (a^) 2 ) 
i(0 2 -(a^) 2 -(a^) 2 + r- 2 (a^) 2 ) 



(48) 



The transverse pressure is P T = (T xx + T vv )/2 and the longitudinal pressure is P L = t 2 T w . All 
these quantities, written here for a single classical field configuration, must of course be averaged 
over the Gaussian ensemble of initial conditions derived in section 2. 

When computing the energy momentum tensor, it is important to do so on a lattice that has 
a fine enough spacing in the rapidity direction, as explained in the previous section. The artifacts 
one obtains at large time if this is not the case are illustrated and discussed in the appendix A. 
In the rest of this section, all the numerical results have been obtained on a 40 x 40 x 320 lattice. 

It is also important to subtract the contribution of the vacuum. Indeed, our resummation 
leads to severe ultraviolet divergences in the energy-momentum tensor that behave generically 
like the fourth power of the ultraviolet cutoff (here, the inverse lattice spacing). This unwanted 
term is a pure vacuum contribution and it can be removed by subtracting the result from another 
computation in which the background field ip in eq. (27) is set to zero ("vacuum"). In all the 
quantities presented later on, this subtraction has been performed. 

5.2 Equation of state 

The first thing to assess is whether the system obeys an equation of state. Towards that end, 
we plot in the figure 6 the energy density and the sum of the three pressures. As in a system 
with a fixed volume (see [67]), the pressure oscillates rapidly initially, with the oscillations being 
damped and disappearing eventually. When this happens, one sees that 



within statistical errors. 

If one recalls the conservation equation 16 that drives the time evolution of the energy density, 



it is also instructive to compare the time dependence of the energy density with two power laws 
that have a special physical meaning: 

i. r^ 1 , expected for a system whose longitudinal pressure is negligible, 

ii. r~ 4 / 3 , expected when P L = P T = e/3. 

As one can see, at the beginning of the evolution, the energy density decreases approximately 
like t -1 , while it is well fitted by r -4 / 3 at later times. This is suggestive of the fact that the 
system is initially highly anisotropic, and then becomes nearly isotropic later on. 

16 This is an exact equation, valid whether the system is in equilibrium or not. Note that since the energy 
density has a smooth derivative at r = 0+ (see for instance the figure 2), the longitudinal pressure must be 
negative and equal to P L = — e at r = + . This is indeed the case but not obvious from the figures, where we 
plot the absolute value of the pressures. 



e = 2P T + P, 



L 1 



(49) 
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Figure 6: Time evolution of the energy density and of the trace of the pressure tensor, compared 
to the power law behaviors r -4 / 3 and r _1 . Note: we are plotting the absolute value of the 
pressure, since at early times its sign changes periodically. 



5.3 Isotropization 

Given the above results on the time dependence of the energy density, let us now look at the 
evolution of the transverse and longitudinal pressures separately. This is illustrated in hgure 
7. (We also plot the energy density and the sum of the three pressures in order to compare 
the timescales for the relaxation of the pressures and for isotropization.) In this plot, one can 
distinguish three stages in the time evolution of the pressures: 

i. At early times (on timescales comparable to that of the relaxation of the pressure or 
shorter), one sees that the magnitude of the longitudinal pressure drops very quickly. 
When well defined quasiparticles can be identified in the spectrum of the theory, this can 
be understood from the kinetic theory formula for the longitudinal pressure 

p ^Im^wi m - (51) 

as a consequence of the red-shifting of the longitudinal momenta due to the expansion 
of the system. Indeed, particles with a non-zero p z eventually escape from the unit slice 
of rapidity whose evolution we are considering. This is illustrated in figure 8. At early 
times, the system we are considering is arguably not amenable to a description in terms 
of quasiparticles. However, the above argument about redshifting applies equally to the 
longitudinal pressure of a system of fields, and implies a rapid decrease of its magnitude. 

ii. Next, a dramatic change of behavior occurs. The longitudinal pressure increases rapidly 
and approaches the transverse pressure. This change begins when the oscillations of the 
pressure have become small. Moreover, from figure 5, we see that this occurs when the 
particle distribution starts to expand and occupy modes in the v direction. 
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50 100 150 200 250 300 

T 

Figure 7: Time evolution of the diagonal components of T piy , and of the trace of the pressure 
tensor. Note: we are plotting the absolute value of the pressures, since at early times their sign 
changes periodically. 




Figure 8: Evolution of the distribution of longitudinal momentum in a slice of rapidity, for 
particles undergoing free streaming. The thick arrows represent the velocities of the particles. 
The thin straight lines represent the trajectories of free particles-at constant velocity. At the time 
n , we assume that the rapidity slice under consideration contains particles with all longitudinal 
momenta. After some period of free streaming to the time T2, only particles whose momentum 
rapidity y equals the space-time rapidity rj are left. 
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iii. Eventually, the longitudinal pressure becomes very close to the transverse one, albeit at 
times much later than the relaxation of the pressure. This confirms what was guessed on 
the basis of the time dependence of the energy density alone. 




Figure 9: Time evolution of the transverse and longitudinal gradients that enter in the pressures. 
The dotted line means that the corresponding quantity is negative. (We have plotted here its 
absolute value.) We also plot the time evolution of the transverse and longitudinal pressures for 
easy reference. 

The difference between the transverse and longitudinal pressures essentially arises from the 
sign with which the various spatial gradients of the field enter in the formulas (48). In particular, 
it easy to see from these formulas that the two pressures are equal provided that 

\(V^) 2 = ^(d v <j>) 2 ■ (52) 

In order to better understand the evolution of the two pressures, we have represented in figure 
9 the expectation values of the two quantities (Vj_</>) 2 /2 and t~ 2 (d v (f>) 2 . Note that these quan- 
tities have been vacuum-subtracted, as explained after eq. (48), which means that they are not 
necessarily positive. Negative values are indicated in the plot by a dotted line instead of a solid 
one. 

At very early times, the transverse gradient is moderately small and positive, while the 
longitudinal gradient is large and negative (after we have subtracted the contribution from the 
vacuum). Up to r « 10, the transverse gradient remains roughly constant (more precisely, it 
oscillates around a constant value), while the magnitude of the longitudinal gradient decreases 
very rapidly. (Note that since it is negative, this means that it is in fact increasing with time.) 

For a short period of time after r rs 10, the transverse gradient increases mildly. From the 
plots showing the evolution of the occupation number in the figure 5, this is clearly related to 
the expansion of the particle distribution in fc j_ . During this period, the transverse gradient is 
much larger than the longitudinal one (which is still negative), and therefore the longitudinal 
pressure is much smaller than the transverse one. 
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A qualitative change in behavior occurs around r ~ 50: the longitudinal gradient evolves 
faster than the transverse one and it increases rapidly, approaching zero to become positive 
shortly after r ss 100. The consequence on the longitudinal pressure is a rapid increase despite the 
expansion of the system. Microscopically the increase of the longitudinal gradients is presumably 
due to the resonance that can reshuffle efficiently momenta now that the expansion rate of the 
system is lower. (Before that point, any particle that was scattered at a non-zero p z immediately 
escaped from the rapidity slice under consideration.) This trend continues until the transverse 
and longitudinal gradients become comparable, at which point the two pressures are nearly equal. 

In appendix B, we study how the isotropization timescale varies when we change the amplitude 
of the background field and the coupling strength; our results show that isotropization is faster 
for larger fields and/or larger couplings. 



6 Comparison with hydrodynamics 

The previous results indicate that the trace of the pressure tensor relaxes to its equilibrium value. 
Subsequently, the longitudinal pressure approaches the transverse one and the system becomes 
isotropic despite the continuing longitudinal expansion. A natural question to ask is whether the 
details of this time evolution are compatible with hydrodynamics 17 . Of course it does not make 
sense to try answering this question before the time at which the pressures become well defined, 
namely, before the pressure oscillations have disappeared. Further, it is clear that to map our 
pressures onto hydrodynamics, a minimal requirement is to include a shear viscosity that will 
allow the pressure tensor to be anisotropic. Given the geometry of our setup, it is sufficient to 
consider boost invariant hydrodynamics. Also, since our problem is translationally invariant in 
the transverse directions, we can set the transverse fluid velocity to zero. 



6.1 Energy and momentum conservation 

In this situation, the time evolution of the energy density is related to the longitudinal pressure, 

A first consistency check that we can make is to verify that this equation holds from our numerical 
results. This is shown in fig. 10, where one clearly sees that this equation is satisfied. This should 
of course not be a surprise because the equation comes directly from the conservation of energy 
and momentum and should be satisfied regardless of whether hydrodynamics applies or not. 



6.2 Comparison with first order viscous hydrodynamics 

The previous check was merely a verification that our numerical approach fulfills basic conserva- 
tion laws. In order to go further, let us recall that in first order viscous Bjorken hydrodynamics 
the longitudinal and transverse pressures would be given by 

e 2n 

P = - H - 

3 + 3r 



17 See [76] for a recent review on relativistic hydrodynamics. One can also find an interesting comparison 
between an exactly solvable (via the AdS/CFT correspondence) toy model at strong coupling and relativistic 
hydrodynamics in ref. [77]. 
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Figure 10: Comparison between the two sides of eq. (53). 



where rj is the shear viscosity. This ansatz assumes that we are already in the regime where 
e = 2P T + P L , and attributes the difference between the two pressures to the effect of the shear. 

First of all, we can compare the time evolution of the pressure anisotropy, (P T — P L )/e, shown 
in fig. 11, with what one would expect in a hydrodynamical expansion. From the time where 
the oscillations of the pressures have died out to the time where this anisotropy has become very 
small, it can be fitted by an exponential form exp(— r 2 /2a 2 ). At later times, the quality of the fit 
can be improved slightly by adding a term in r -2 / 3 . (Note that at times greater than r > 280, 
the anisotropy rises again due to lattice artifacts.) 

If the decrease of the anisotropy was driven by viscous hydrodynamics, one would expect 
its behavior to be a power law r -2 / 3 , not an exponential. Although it is not easy to pinpoint 
what precisely is causing this exponential falloff of the anisotropy, it seems very plausible that 
the presence of instabilities in the microscopic dynamics of the system is responsible. One may 
argue that the hydrodynamical regime only starts when the power law term in r~ 2 / 3 becomes 
important -roughly after r ~ 200 lattice spacings. Then one could turn the coefficient b of the 
r~ 2 / 3 term into a value of the ratio rj/s. The dimensionless ratio rj/s determines the strength 
of the viscous effects on the hydrodynamical evolution of the system. In a scale invariant theory 
like the one we consider here, it is natural that this ratio is a fixed number 18 that depends only 
on the value of the coupling constant g. In the hydrodynamical regime, one would have 

f ~ S Tf ~ s 41/4 T 2 /3 ' 1 ' 

e J hydro S Te >° ™ ' j T ' 

b 

where we have used the Stefan-Boltzmann formulas to estimate the energy and entropy density 
(assuming the system is close to equilibrium), 

^ 2?r2T3 ~ 3/4 

£= ^Q- ' S= ^5~ ' S " £ ' (56) 

18 This statement is true at the order at which we analyze the system but would be violated once one includes 
running coupling effects. 
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Figure 11: Time evolution of the pressure anisotropy (P T — P L )jt. Red dotted line : fit in 
exp(— r 2 /2a 2 ). Black solid line : fit in exp(— t 2 /2a 2 ) plus a power law. The relaxation time in 
the exponential is a 82.9 lattice spacings, and the coefficient of the power law correction is 
b ss 0.182 (in lattice units). 



and where A is the coefficient in the asymptotic behavior of the energy density, e 
From this formula, we extract the following value for the ratio r)/s, 



0.26 



At- 4 / 3 . 



(57) 



with a large systematic uncertainty due to the limitations of our analysis outlined previously. 

To go beyond this simple comparison, we should solve eq. (53) with the ansatz of eq. (54) for 
the longitudinal pressure. To close the equation, we need a way to relate rj to e, with eq. (56). 
Even without solving the hydrodynamical equations, we can use this approximation for the 
entropy density in order to extract an effective rj/s from our numerical results. This amounts to 
rewriting eqs. (54) as 



2 

37 



,3/4 



err 



p =i-A 

L 3 3r 



,3/4 



c rr 



(58) 



Contrary to the estimate in eq. (57), where we attributed to the viscosity only the power law 
term in the difference P T — P L , we now define this effective viscosity from the full difference, 
including also the exponential term. This effective ratio can be computed at each time, and is 
shown in the figure 12. The first thing one sees is that it is not constant, which certainly means 
that there are discrepancies with hydrodynamics. The value of this ratio before r ~ 70 probably 
does not make much sense, because the pressures are oscillating before this time. After that, the 
ratio starts near a value of 10 and decreases rapidly to reach a value of about 0.5. In addition 
to the non-constancy of this ratio, it should be noted that its value is quite small compared to 
the perturbative value. For a real scalar field theory with a g 2 (f> 4 /4\ coupling, the lowest order 
perturbative result [78] is roughly 

5-1$. (59) 
s g 4 
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Figure 12: Evolution of the effective rj/s ratio. We have also represented the perturbative value 
of this ratio, the value extracted from the power term (eq. (57)), and the conjectured lower bound 
1/47T derived in AdS/CFT at strong coupling. 



Given that g — 4 in our simulation, one would expect r\j s ~ 40, which is much larger than 
the effective ratio we observe. In other words, the system of fields we have studied is a much 
better fluid than one would expect on the basis of the perturbative estimates of the viscosity 
to entropy density ratio. This anomalously small value of the r\j s ratio is perhaps related to 
the phenomenon discussed in [79, 80] where it is argued that, in systems subject to turbulent 
unstable fields, momentum transport may occur as if the viscosity were much smaller than naive 
estimates from transport cross-sections. 

To carry the comparison of our simulations with hydrodynamics even further, we must solve 
the hydrodynamic equations; using eq. (56), the evolution of e in viscous hydrodynamics has the 
closed form 

OT 3 t 3s T 

Comparisons between the two frameworks can be performed as follows: 

i. Choose an initial time To, where hydrodynamics is initiated. This time should be well after 
the oscillations in the pressures have disappeared. 

ii. Initialize the energy density so that it has the same value in the two frameworks at To- 

iii. Set the value of the ratio r]/s in order to have the right values for the transverse and 
longitudinal pressures at t . This is always possible if To is in the region where e = 2P T +P L ■ 

iv. Solve (numerically) the differential equation (60) from these initial conditions and compare 
with the evolution one obtains in classical statistical field theory. 

The result of this comparison is shown in fig. 13, for three values of the initial time To for 
hydrodynamics. The main lesson from this figure is that the pressure tensor isotropizes much 
faster in classical statistical field theory than it does in hydrodynamics. Hydrodynamics begun 
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Figure 13: Comparison between classical statistical field theory and first order viscous hydrody- 
namics, for various starting times of the hydrodynamical evolution (tq = 70, 100, 200). Solid lines: 
pressures in classical statistical field theory. Dots: results of first order viscous hydrodynamics. 

at earlier times cannot compete and thus flow and viscosity bounds derived from such exercises 
are suspect. 

The last pair of hydro curves, with an initial time tq — 200, remains close to the true 
evolution, but this is only because it starts with an already small amount of anisotropy. Therefore, 
for hydrodynamics to describe accurately the evolution, hydrodynamical evolution should be 
initiated at a time where the pressure tensor is already close to isotropy. Our comparison is 
based on first order viscous hydrodynamics only. Arguably, viscous corrections are large in the 
system when P T and P L differ substantially and it is unclear whether first order hydrodynamics 
is a valid approximation. 

We close this section with a generic comment about the isotropization time. The precise 
value of the time we obtain in this model of scalar fields is of little relevance for the analogous 
question in QCD because the two theories are quite different in many respects. However, we 
expect the conclusion that the true isotropization time is considerably shorter than the time one 
would infer from hydrodynamics (or equivalently from a Boltzmann equation) is quite generic 
in theories with instabilities. The other lesson from this comparison is that such systems may 
have an effective viscosity to entropy density ratio which is much smaller than what one expects 
from perturbation theory. Therefore one could have nearly ideal fluid behavior without invoking 
strong coupling. Moreover, we know from [68] that full thermalization of this system takes much 
longer than the times considered here. We may conclude that full thermalization is not necessary 
for a hydrodynamical description; only near isotropy is important. 

7 Summary and outlook 

This work extends our previous study of the role of quantum fluctuations on the thermalization 
of strong fields (subject to instabilities) to the situation where the system undergoes longitudinal 
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expansion. We considered a one-component scalar field theory, which is much simpler to study 
numerically than Yang-Mills theory, and implemented the same resummation program that one 
would employ in studying the early stages of heavy ion collisions in the Color Glass Condensate 
framework. 

As in the case of a system confined to a fixed box, the trace of the pressure tensor exhibits 
rapid temporal oscillations that are quickly damped via the phenomenon of phase decoherence. 
Following this state we find that the sum of the three pressures become equal to the energy 
density (a well-define equation of state exists), but the transverse and longitudinal pressures 
remain different until full isotropization at an even later time. 

The main new result of this study is that the pressure tensor eventually becomes isotropic 
despite the longitudinal expansion of the system. After an initial decrease due to the rapid 
expansion at early times, the longitudinal pressure increases at later times to reach a value that 
is very close to that of the transverse pressure. This goes a long way towards justifying the 
validity of the hydrodynamical description for expanding systems, as encountered in heavy ion 
collisions. 

We also studied the time evolution of the particle distribution. Even with a very far from 
equilibrium initial condition that has a single k± mode at t = + , a continuum of modes 
rapidly get filled by the nonlinear interactions of the fields. This study of the occupation number 
highlights an important limitation of numerical simulations that use rapidity as the longitudinal 
coordinate, namely, the particle distribution expands in v -the conjugate momentum of rapidity- 
over time and eventually reaches the ultraviolet cutoff imposed on this variable by the lattice 
spacing. Therefore in order to pursue the study until sufficiently large times without being 
affected by these lattice artifacts, it is necessary to use a lattice that has a very fine mesh in the 
rapidity direction. This issue is by no means specific to a scalar theory, and will be present when 
this approach is applied to Yang-Mills theory. 

In all the plots presented in this paper, all dimensionful quantities were expressed in units 
of the appropriate power of the transverse lattice spacing. Since the </> 4 scalar theory in 3 + 1 
dimensions is scale invariant at the classical level (as is Yang-Mills theory), one can re-express 
all of them in terms of a single dimensionful physical parameter. In the Color Glass Condensate 
framework, this parameter would be the saturation momentum Q s of the colliding nuclei. Doing 
this for QCD would provide an answer to the following fundamental question: how long does 
isotropization take, given the value of the saturation momentum and of the strong coupling 
constant? However, the lessons of this exercise in scalar field theory cannot be simply applied to 
the gauge theory case because the former differs from QCD in a number of crucial ways. To list 
a few, (i) the field content of QCD is much richer, (ii) there is no meaningful way to compare the 
coupling constants of the two theories, (iii) the strength of their instabilities arc quite different, 
both in their range and in their growth rate. 

Thus it is best to see the present work as a proof of the concept that instabilities can isotropize 
the pressures of an expanding system. The next step obviously is to apply the same treatment to 
Yang-Mills theory. In a previous work [66] , we derived the spectrum of initial field fluctuations 
that one must use in this computation. The numerical implementation of this approach to a 
Yang-Mills theory is under way but will require a significantly larger computational effort than 
what was needed for the present scalar theory. 
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A Lattice artifacts at late times 

As we have already explained in section 4.2, the time extent of the simulation is limited by the 
finite size of the lattice. We demonstrated that as the time increases the particle distribution 
expands to larger v values, [y is the momentum conjugate to rapidity), until it is artificially 
cut-off by the discretization in rapidity. 




P L /e: N = 80 160 320 — 

P x /e: N = 80 160 320 



Figure 14: Time evolution of the ratios of longitudinal and transverse pressure to the energy 
density, for several choices of discretization in rapidity. 

In figure 14, we show the practical significance of this effect on the numerical computation of 
the components of the energy-momentum tensor, for lattices with increasing number of longitu- 
dinal spacings N (all describing a unit slice of rapidity) . We see that when the number of lattice 
spacings in the rapidity direction is small (see the curves for the 40 x 40 x 80 lattice), the ratio 
P L /e approaches the isotropic value 1/3 and then departs from this value. Likewise, the ratio 
P T /e approaches 1/3 from above, before deviating from this value (in such a way that 2P T + P L 
remains equal to e). The time at which this occurs agrees with the time when the particle distri- 
bution reaches the longitudinal momentum cutoff as shown in table 1. Figure 14 demonstrates 
that this artifact appears at later times as the number of lattice spacings in rapidity is increased. 

The additional complication of renormalization must also be considered when simulating 
scalar fields on a lattice. While there is no mass term in our bare Lagrangian, vacuum fluctuations 
induce a non-zero mass. In a scalar </> theory, the leading contribution stems from the tadpole 
graph having a quadratic UV divergence. The value of this tadpole is given as 
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Figure 15: Tadpole contribution to the mass as a function of time, for various longitudinal lattice 
spacings (and a coupling constant g = 4). 



which clearly exhibits the nature of its divergence. On the lattice the above integral is replaced by 
a discrete sum having a UV cutoff dictated by the size of the lattice spacing. In t-t\ coordinates 
this mass is r dependent (a fixed UV cutoff in v corresponds to a time dependent cutoff in 
k z ~ v l T \ Figure 15 shows the time evolution of m 2 for various discretizations in rapidity. At 
intermediate times the mass has a logarithmic sensitivity to the lattice spacing. 

In practice this means that computations done with the same bare Lagrangian but having 
different longitudinal discretizations correspond to different renormalized theories. The loga- 
rithmic dependence of this mass with N means that it varies more rapidly at small N than at 
large AT, something that can be seen in the figure 14 by comparing the changes from N = 80 
to N = 160 and from N — 160 to N = 320. In the present work, we have not attempted 
to renormalize the bare parameters as this will be further complicated by the fact that we are 
treating the transverse and longitudinal coordinates on different footings. Therefore, one should 
keep in mind that the theory for which results are presented in the paper is not a truly massless 
theory, and the comparison between computations done at varying lattice spacings can only be 
made qualitatively. We should stress that this issue should not arise in gauge theories. In this 
case, the gauge invariance of the Wilson action prevents the generation of a mass term for the 
gluons. 

B How generic is it? Varying g 2 and ip 

A natural question that arises is whether the isotropization that we have observed is a generic 
phenomenon that occurs for all choices of coupling strengths and background fields. 

The scale invariance (at the classical level) of the theory we are studying suggests that it 
should be generic for all background fields since a larger background field can be hidden in a 
rescaling of the time and spatial coordinates. Due to the limitations inherent to the lattice setup, 
we can only vary the background field in a limited window; nevertheless, we have indeed observed 



27 



50 100 150 200 250 300 

X 

Figure 16: Dependence of the isotropization of the pressure on the amplitude of the background 
field. The coupling strength is the same in the four computations. 

numerically that isotropization occurs regardless of the amplitude of the background field. This 
is shown in the figure 16, where the time dependence of (P T — P L )/e is represented for uniform 
background fields of amplitudes ifo = 15, 30, 50, 75. Moreover, we see that the isotropization time 
decreases if <po increases. This decrease is expected because all time scales should be inversely 
proportional to the field amplitudes. Note that there are important lattice artifacts for the 
smallest value of ifo that are visible in the very small change of the isotropization curve when 
going from (p = 15 to ip = 30. In this regime of smaller fields, the tadpole mass-due to the 
vacuum fluctuations that we have not renormalized-seems to dominate the dynamics. This is a 
lattice issue related to the breaking of the scale invariance by lattice cutoffs of the theory we are 
simulating. 

Less obvious is the question of what happens when varying the coupling. We do expect the 
isotropization time to decrease if we increase the coupling; however, there could be some critical 
coupling below which isotropization never occurs. What we observe numerically over the range 
of couplings explored (limited by the lattice configurations) is that isotropization is faster at 
larger coupling strengths, as shown in fig. 17. 

C Computation of Hankel functions of imaginary index 

In the computation of the initial conditions (eq. (27)) and of the occupation number (eq. (45)), 

(2) 

we need Hankel functions of imaginary index H iv . They are both needed at small arguments 
(for the initial conditions), and at fairly large arguments (for the occupation number). And we 
need them for very large indices v, since on a lattice with N spacings in a unit rapidity interval, 
the maximal value of v is ^ max = 2N '. Thus, for the value N — 320 that we have used in most 
of this paper, we need Hankel functions with indices up to 640. 

To the best of our knowledge, Hankel functions of imaginary index are not implemented in 
the common scientific numerical libraries such as CERNLIB or the GNU Scientific Library, but 
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Figure 17: Dependence of the isotropization of the pressure on the strength of the coupling. 
The amplitude of the background field is the same in the three computations. 

it is fairly easy to evaluate them numerically with high accuracy. First of all, one should recall 
that they are solutions of the Bessel equation 

/" + -/' + (1 + 4)/ = °- (62) 

Naturally, a pair of linearly independent solutions of this equation can be constructed numerically 
to any desired accuracy. But the problem is to find the proper linear combination of those that 
gives the desired Hankel function. It is best defined by its asymptotic expansion, 




This expansion should be used with caution, because it is not a convergent series if summed to 
arbitrarily large n. But at a given v 1 one can find a sufficiently large r and an optimal n such 
that the residual term is extremely small. In practice, one computes the successive terms in 
eq. (63) for a given v and a large r: they decrease for small k, but eventually start increasing 
again for k larger than some n. The last computed term before they start increasing is used as 
an estimate of the error. If it is satisfactorily small, then we compute the Hankel function and 
its first derivative by summing the asymptotic expansion up to this n. If the error is not small 
enough, we start over this process for the same v, but a larger t. 

At this point, we know the value of H iv ' (and that of its first derivative) for some value of r. 
The next step is to use these as initial conditions to solve numerically the Bessel equation, both 
forward and backward in order to cover all the desired range of arguments. This can be done to 
very high accuracy by using a high order solver with adaptive steps. 
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